##################################################################
# This file runs the main and supplementary analyses.
##################################################################
#set working directory
setwd("/Users/agunderson/Dropbox/Gender and Policing/Female Cops Paper/PoP")
#clear working environment
rm(list=ls())
#load packages
library(ggplot2)
library(tidyverse)
library(lfe)
library(stargazer)
library(xtable)
library(naniar)
library(readxl)
library(cdlTools)
options(scipen=999) #to allow merge in of very large numbers

##################################################################
# Model functions
##################################################################
main_model_html <- function(dataframe,tabletitle,label,dv1,dv2,iv1,ivlab,depvarlabel){
  model1 <- felm(dv1~iv1
                 |year|0|0,data=dataframe)
  
  #model 2: female police officers and controls
  model1controls <- felm(dv1~iv1+totalemployees+log(budget)+percblackfulltime+raperate+sheriff+communitypolicing
                         |year|0|0,data=dataframe) 
  
  #Model 2: female police leadership with diff controls 
  model2controls <- felm(dv1~iv1+totalemployees+ percblackfulltime+raperate+sheriff
                         |year|0|0,data=dataframe)
  #with reports
  model1a <- felm(dv2~iv1|year|0|0,data=dataframe)
  
  #model 2: female police officers and controls
  model1acontrols <- felm(dv2~iv1+totalemployees+log(budget)+percblackfulltime+sheriff+communitypolicing
                          |year|0|0,data=dataframe) 
  
  #Model 2: female police leadership with diff controls 
  model2acontrols <- felm(dv2~iv1+totalemployees+percblackfulltime+sheriff|year|0|0,data=dataframe)
  
  #Export for Latex -----
  stargazer(model1a,model1,model2acontrols,model2controls,model1acontrols,model1controls,covariate.labels = c(ivlab,
                                                                                                              "Total Agency Employees","Logged Agency Budget",
                                                                                                              "Perc. Black Officers","Rape Report Rate","Sheriff",
                                                                                                              "Community Policing"),
            label=label,
            title=tabletitle,
            dep.var.labels=depvarlabel,
            style="apsr",notes="Year dummy variables included.",
            out= paste0(label,".html"),
            font.size = "tiny",table.placement="h!",column.sep.width = "-10pt")
}

##################################################################
# Import data --------------
################################################################## 
genderpolice <- read_csv("PoPData.csv")%>%
  mutate(rapeclearancerate = (rapescleared*100000)/population) %>%
  na_if(., Inf) %>%
  group_by(ORI) %>%
  filter(n()<10) %>% #some mismatches in names, so need to exclude those from analysis (more than one agency matched on name only)
  filter(percblackfulltime<101)%>% #two agencies have more black employees than total employees
  as.data.frame() 

##################################################################
# Figures
################################################################## 
ggplot(genderpolice, aes(x=factor(year), y=percfemale)) + 
  geom_boxplot()+theme_bw()+ylab("Percent Female Sworn Officers")+
  xlab("Year")
ggsave("Figure1.pdf") 

plot <- genderpolice %>% dplyr::select(percfemale,sumrapearrestsrate,raperate) %>%
  gather(., Variable, Rate, sumrapearrestsrate:raperate, factor_key=TRUE) %>% 
  mutate(Variable = str_replace_all(Variable, c("sumrapearrestsrate" = "Rape Arrest Rate", "raperate" = "Rape Report Rate"))) 

ggplot(plot,aes(x=percfemale,y=Rate,shape=Variable,linetype=Variable,color=Variable))+geom_point(alpha = 0.3)+stat_smooth(method="lm",se=F)+
  scale_y_continuous(trans="log10")+theme_bw()+ scale_colour_grey(start = 0, end = .75) + xlab("Percent Female Officers")
ggsave("Figure2.pdf") 

#################################################################
# Tables
##################################################################
#table 1
stargazer(genderpolice[c("percfemale","femalefulltime","sumrapearrestsrate","totalemployees","percblackfulltime","raperate","sheriff",
                         "budget","communitypolicing")],
          covariate.labels = c("Perc. Female Officers","No. Female Officers", "Rape Arrest Rate","Total Agency Employees",
                               "Perc. Black Officers", "Rape Report Rate","Sheriff's Office","Agency Budget","Community Policing"),
          label="summarystats",digits=2,
          title="Gender Composition and Crime Data of U.S. Police Forces, 1987-2016",
          out="Table1.html",
          font.size = "scriptsize",table.placement="h!")

genderpolice <- genderpolice %>% filter(sumrapearrestsrate<1000) #filter out massive outlier

#table 2
main_model_html(genderpolice,"The Effect of Female Police on Rape Report and Arrest Rates, 1987-2016",
                "Table2",genderpolice$sumrapearrestsrate,genderpolice$raperate,genderpolice$percfemale,
                "Perc. Female Officers",
                c("Rape Report Rate","Rape Arrest Rate","Rape Report Rate",
                  "Rape Arrest Rate","Rape Report Rate","Rape Arrest Rate"))
#table 3: violent crime
model1 <- felm(sumviolarrestsrate~percfemale
               |year|0|0,data=genderpolice)
model1controls <- felm(sumviolarrestsrate~percfemale+totalemployees+log(budget)+percblackfulltime+violentcrimerate+sheriff+communitypolicing
                       |year|0|0,data=genderpolice) 
model2controls <- felm(sumviolarrestsrate~percfemale+totalemployees+ percblackfulltime+violentcrimerate+sheriff
                       |year|0|0,data=genderpolice)
model1a <- felm(violentcrimerate~percfemale|year|0|0,data=genderpolice)
model1acontrols <- felm(violentcrimerate~percfemale+totalemployees+log(budget)+percblackfulltime+sheriff+communitypolicing
                        |year|0|0,data=genderpolice) 
model2acontrols <- felm(violentcrimerate~percfemale+totalemployees+percblackfulltime+sheriff|year|0|0,data=genderpolice)

stargazer(model1a,model1,model2acontrols,model2controls,model1acontrols,model1controls,covariate.labels = c("Perc. Female Officers",
                                                                                                            "Total Agency Employees","Logged Agency Budget",
                                                                                                            "Perc. Black Officers","Violent Crime Report Rate","Sheriff",
                                                                                                            "Community Policing"),
          label="femalecopsviolentcrimerate",
          title="The Effect of Female Police on Violent Crime Arrest and Report Rates, 1987-2016",
          dep.var.labels=c("VC Report Rate","VC Arrest Rate","VC Report Rate",
                           "VC Arrest Rate","VC Report Rate","VC Arrest Rate"),
          style="apsr",notes="Year dummy variables included.",
          out= "Table3.html",
          font.size = "tiny",table.placement="h!",column.sep.width = "-10pt")
